Modulate stress distribution with bio-inspired irregular architected materials towards optimal tissue support

Natural materials typically exhibit irregular and non-periodic architectures, endowing them with compelling functionalities such as body protection, camouflage, and mechanical stress modulation. Among these functionalities, mechanical stress modulation is crucial for homeostasis regulation and tissue remodeling. Here, we uncover the relationship between stress modulation functionality and the irregularity of bio-inspired architected materials by a generative computational framework. This framework optimizes the spatial distribution of a limited set of basic building blocks and uses these blocks to assemble irregular materials with heterogeneous, disordered microstructures. Despite being irregular and non-periodic, the assembled materials display spatially varying properties that precisely modulate stress distribution towards target values in various control regions and load cases, echoing the robust stress modulation capability of natural materials. The performance of the generated irregular architected materials is experimentally validated with 3D printed physical samples — a good agreement with target stress distribution is observed. Owing to its capability to redirect loads while keeping a proper amount of stress to stimulate bone repair, we demonstrate the potential application of the stress-programmable architected materials as support in orthopedic femur restoration.

Supplementary Figure 1: Material database creation in two-and three-dimensional cases.This creation process involves specifying basic building blocks, sampling frequency combinations, generating disordered microstructures, and evaluating microstructural homogenized elastic moduli.

Prescribing basic building blocks
Basic building blocks are fundamental ingredients for creating disordered microstructures, and these microstructures further compose the irregular architected materials.To ensure the diversity and seamlessness of the generated materials, we require that the basic building blocks are representative and can be connected via rotation.Specifically, we prescribe four types of basic building blocks, namely cross, arrow, corner, and line shapes, in the two-dimensional (2D) case and four types of basic building blocks, namely cross, corner, t, and cross-line shapes, in the three-dimensional (3D) case.These specified building blocks exhibit diverse geometries with each block capable of connecting only in specific directions to others, achieving a large material property space (see Supplementary Figure 3 that will be introduced later).

Sampling frequency combinations of microstructures
We denote one frequency combination of a disordered microstructure by {ξ 1 , ξ 2 , ξ 3 , ξ 4 } satisfying 4 i=1 ξ i = 1, where ξ 1 -ξ 4 are the frequencies of the four basic building blocks in the microstructure, respectively.To create a representative material database, we uniformly sample 200 microstructural frequency combinations on the hyperplane 4  i=1 ξ i = 1 for both 2D and 3D cases following the procedure in the Supporting Information of [1].

Generating disordered microstructures
Based on the prescribed basic building blocks and sampled frequency combinations, we use the original virtual growth scheme [1] to create 100 specimens of disordered microstructures for each frequency combination, resulting in 20,000 square specimens in 2D and 20,000 cubic specimens in 3D.Each 2D specimen contains 40 × 40 building blocks, and each 3D specimen contains 10 × 10 × 10 building blocks.

Evaluating microstructural homogenized elastic modulus
The microstructure generated by the virtual growth scheme is irregular and should be homogenized to evaluate its elastic modulus.We use a numerical homogenization approach [2] (see also Section 1.2 for details) to evaluate the microstructural homogenized elastic modulus -a 3 × 3 matrix in 2D (D (2) ) and a 6 × 6 matrix in 3D (D (3) ), expressed in matrix notation [3].We then compute the average elastic modulus for every 100 specimens corresponding to the same frequency combination.Finally, we obtain 200 pairs of frequency combination-elastic modulus relationships for both 2D and 3D cases.

Numerical homogenization of disordered microstructures
During the material database creation, we need to evaluate the homogenized elastic modulus (D) of disordered microstructures.We achieve this goal by leveraging a numerical homogenization approach in [2].This approach assumes that the microstructure under consideration is a unit cell of an infinitely large periodic structure, and we can compute the homogenized elastic modulus by imposing periodic boundary conditions on the unit cell (disordered microstructure herein).For the reader's convenience, we outline the key steps of this approach as follows.
We begin by specifying the periodic vectors, a 1 , a 2 , and a 3 , of the irregular microstructure.Next, we discretize the microstructure into finite elements and classify the degrees of freedom (DOFs) as independent or dependent.We then find the topological matrices, B 0 and B a , such that where d and d 0 are the displacement vectors of all DOFs and independent DOFs, respectively.The variable Δa = [Δa 1 , Δa 2 , Δa 3 ] contains the changes in periodic vectors due to deformation, namely Δa 1 , Δa 2 , and Δa 3 .Note that Eq. (1) essentially represents the periodic boundary condition, and it implies a separation of lengths scales between the unit cell and the material piece.
Based on the finite element analysis (FEA), we evaluate the global stiffness matrix of the discretized microstructure as K uc .In this study, we use first-order quadrilateral elements in 2D and frame elements in 3D to improve computational efficiency.By applying the periodic boundary conditions, we can establish the periodic equilibrium equation as Substituting (1) into (2) yields and we obtain the displacement vector of independent DOFs as Substituting ( 3) into (1) gives the displacement vector of all DOFs as We then evaluate the stored energy density as We apply a dummy macroscopic strain field of ε M,3×3 on the microstructure.According to the definition of the infinitesimal strain, we have Δa . Finally, we can rewrite the stored energy density as Based on continuum mechanics [3], the homogenized elastic modulus (in matrix notation) is given as and we recall , where B 0 , B a , B ε , K uc , and V are ready to compute for a given microstructure.

Material property distribution
Considering the stochastic nature of the virtual growth algorithm, the material properties in the generated database exhibit randomness.In this subsection, we study the material property distribution of the created specimens.Supplementary Figure 2 displays the distribution of the D 11 component of the homogenized elastic modulus (D) for 100 specimens corresponding to randomly selected 9 frequency combinations.In Supplementary Figure 2, the x-axis represents the material property (D 11 in this case).The left y-axis denotes the statistical frequency (count), while the right y-axis represents the probability density function.The inset in each panel presents the associated frequency combination of building blocks.Notably, the material properties of the generated specimens generally exhibit unimodality and near-zero skewness.Therefore, we assume that the material properties follow a normal distribution with their mean values determined by the frequency combinations of building blocks.

Achieved material property space
In this subsection, we investigate the attained material property space based on the discrete material database in Supplementary Figure 1.Following the convention in [1], we begin by defining the average Young's modulus and Poisson ratio to quantify the homogenized anisotropic material properties of the microstructures in the database.Specifically, we define the average Young's modulus as and the average Poisson's ratio as where S represents the fourth-order homogenized compliance tensor that can be derived from the matrix D.
With the above definitions, we plot the material property space for both 2D and 3D scenarios in Supplementary Figure 3. Here, the x-axis is the relative average Young's modulus defined as E ave /E solid , and E solid is the Young's modulus of the solid material; the y-axis is the average Poisson's ratio, ν ave .In Supplementary Figure 3, the color wheels signify the frequency combinations present in the material database in Supplementary Figure 1, and their positions correspond to E ave /E solid -ν ave pairs averaged over the specimens of the same frequency combination.The shaded area in Supplementary Figure 3 indicates the actual material properties achieved by all the specimens.The insets in both Supplementary Figures 3a and  b illustrate specimens generated based on their frequency combinations, respectively.According to Supplementary Supplementary Figure 3, the 2D material property space spans E ave /E solid values in [0.001, 0.111] and ν ave values in [-0.309, 0.365]; the 3D material property space spans the counterparts in [0.003, 0.036] and [0.008, 0.144], respectively.We remark again that the broad spectrum of material properties results from the selection of the representative building blocks.

a b
2D material property space 3D material property space Supplementary Figure 3: Material property space in 2D and 3D cases.The x-axis is the relative average Young's modulus (E ave /E solid ), and the y-axis is the average Poisson's ratio (ν ave ).The color wheels represent the frequency combinations in the material database, and their positions represent the corresponding E ave /E solid -ν ave pairs averaged over specimens of the same frequency combination, respectively.The shaded regions indicate the actual material properties achieved by all the specimens, and the insets illustrate the generated microstructures.

Machine learning model for predicting microstructural elastic modulus
The above material database stores discrete relationships between the frequency combinations and the homogenized elastic modulus of disordered microstructures.However, we need continuous and differentiable relationships for the subsequent topology optimization.To capture such relationships, we train a machine learning model -a fully connected neural network -and use it to predict homogenized elastic modulus based on the input frequency combination.The architecture of the neural network as well as the training and testing processes are outlined below.

Architecture of the neural network
The architecture of the neural network is shown in Supplementary Figure 4.This neural network comprises one input layer with N b = 4 nodes corresponding to the 4 entries of the input frequency combination.It also has two hidden layers with m and p nodes, respectively, and one output layer with q nodes corresponding to the maximum number of independent entries of homogenized elastic modulus (q = 6 in 2D and q = 21 in 3D).We employ the rectified linear activation unit, ReLU(•) = max{•, 0}, as the nonlinear activation function.The mathematical expressions of the architecture in Supplementary Figure 4 are given as

Training and testing the neural network
To optimize the weighting matrices (W 1 , W 2 , and W 3 ) and bias vectors (b 1 , b 2 , and b 3 ) of the neural network, we use the mean squared error (MSE) of the microstructural elastic modulus as the loss function (L) and minimize it using an Adam optimizer [4].To configure this optimizer, we set the coefficients used for computing running averages of gradient and its square as β = (0.900, 0.999) and add a term ε = 10 −8 to the denominator to improve numerical stability.Before training the neural network, we also need to determine its optimal architecture for a faster convergence speed and a smaller value of the loss function, L. To achieve this, we randomly test 100 sets of architecture parameters, including the learning rate of the Adam optimizer, l ∈ {10 −5.0 , 10 −4.5 , 10 −4.0 , 10 −3.5 , 10 −3.0 , 10 −2.5 , 10 −2.0 }, and the numbers of nodes in the two hidden layers, m and p ∈ {4, 8, 16, 32, 64, 128, 256, 512}.We test each set of architecture parameters, (l, m, p), by running the neural network for 200 epochs.We finally select the architecture parameters as l = 10 −2 , m = 512, and p = 256, giving the minimum loss, L, as shown in Supplementary Figure 5a.... Based on the selected architecture parameters of the neural network, we now formally train its weighting matrices (W 1 , W 2 , and W 3 ) and bias vectors (b 1 , b 2 , and b 3 ).Specifically, we randomly set 160 frequency combination-elastic modulus pairs in the previous material database as training data and the remaining 40 pairs as testing data.We train the neural network for 1000 epochs using an adaptive learning rate: 10 −2 for epochs 1-200, 10 −3 for epochs 201-500, and 10 −4 for epochs 501-1000.We implement and train this neural network with an in-house Python code [5] and the PyTorch package [6].After the training, the loss function of the neural network (L) converges to a sufficiently small value (Supplementary Figure 5b).Based on the testing data, the neural network performs well in predicting the microstructural homogenized elastic modulus with small MSEs of 0.134 MPa in 2D and 0.296 MPa in 3D, compared to their corresponding mean values of 2.687 MPa and 7.344 MPa, respectively.

Macroscopic topology optimization for stress modulation
Based on the continuous relationship between the microstructural frequency combination and elastic modulus, we can establish a macroscopic topology optimization formulation for stress modulation.This formulation consists of design space parameterization, elasticity interpolation, and static equilibrium equations, which are introduced below.After that, we introduce the complete formulation and show the detailed problem setups for stress modulation.

Design space parameterization
We need to parameterize the design space to optimize the microstructural layout and frequency combination.Specifically, we parameterize the microstructural layout with a density field, ρ(X) ∈ [0, 1], and parameterize the microstructural frequency combination with frequency fields, where N b represents the number of types of basic building blocks, and ε ∈ [0, 0.5) is a prescribed constant.The variable X is the position vector of a material point in the design domain, Ω.
To reduce the mesh dependency and eliminate the checkerboard pattern of the density field, ρ(X), during topology optimization, we apply a density filter [7] and derive the filtered density field as , where the weighting factor, w(X, X ) = max{0, R − |X − X |}, is determined by the filter radius (R) and the Euclidean distance (|X − X |) between material points X and X .To ensure the 0-1 solution of the density field, we further apply a Heaviside projection [8] and derive the physical density field as where θ = 0.5 is a threshold value and β is a sharpness parameter.The physical density ρ(X) = 1 represents the presence of disordered microstructures, and ρ(X) = 0 represents the voids.
As for the frequency fields, they need to satisfy the axiom of probability, N b i=1 ξ i (X) = 1.To effectively impose this constraint, we propose a normalization projection scheme and obtain the physical frequency fields as where ξ i (X) represents the i-th component of the microstructural frequency combination.

Elasticity interpolation and static equilibrium equations
After parameterizing the design space with the physical density (ρ(X)) and frequency (ξ i (X) for i = 1, 2, ..., N b ) fields, we need to interpolate the microstructural elastic modulus (C(X), in a tensorial form) for these variables.Based on the modified solid isotropic material with penalization (SIMP) [9], the interpolated microstructural elastic modulus is given as where ε ρ = 10 −6 is a small positive number preventing the singularity of the stiffness matrix in FEA.The parameter p ρ = 3 penalizes the intermediate values of the physical density field, ρ ∈ (0, 1), and leads to a 0-1 solution of ρ.The parameter α characterizes the designable region as α = 0 and the passive region as α = 1.The variables C d and C p are the corresponding elastic moduli, respectively, where C d can be predicted by the trained machine learning model according to the inputs of ξ 1 , ξ 2 , ..., ξ N b .
Based on the interpolated microstructural elastic modulus, we establish the static equilibrium equations as [3] ⎧ ⎨ ⎩ Div(C : where the variable E = (∇u + ∇u )/2 is the infinitesimal strain tensor, and u is the displacement field.The variable b is the applied body force, and u is the displacement field prescribed on one part of the boundary, ∂Ω D .The variable t is the traction applied on the remaining part of the boundary, ∂Ω N , with a unit outward normal vector, N.

Macroscopic topology optimization formulation
Based on the design space parameterization and static equilibrium equations, we can employ topology optimization to modulate stress distribution.Topology optimization [10,8] is an approach that determines the material layout that optimizes a given objective function while satisfying certain constraints.Leveraging topology optimization, researchers successfully design optimal structures and materials with complex programmable behaviors [11,12,13,14].In this work, we propose a macroscopic topology optimization formulation for stress modulation, which is expressed as with: static equilibrium equations in ( 4), where J is the objective function representing the stress modulation error accounting for N s stress measures and N l load cases.The parameter w [l] ∈ [0, 1] is the weighing factor of the l-th load case and satisfies is the indicator of stress control regions of the s-th stress measure under the l-th load case, where θ [l,s] = 1 represents stress control regions, and θ [l,s] = 0 represents noncontrolled regions.The variable σ [l,s] (X) is a stress measure, and σ [l,s] (X) is the corresponding target value.The parameter p is used in the p-norm expression to approximate the maximum stress modulation error.The symbol g represents a volume constraint function, and the parameter V represents the allowable volume fraction of disordered microstructures.
To implement this proposed macroscopic topology optimization formulation, we utilize an in-house Python code [5] with the open-source FEniCSx package [15], which features parallel computing and automatic differentiation for the FEA and sensitivity analysis.We also use the method of moving asymptotes (MMA) [16] to update the design variables, ρ and ξ i for i = 1, 2, ..., N b .

Detailed problem setups and additional investigations
Based on the proposed macroscopic topology optimization framework, we move forward to optimize the microstructural layout and frequency combination to minimize the stress modulation error.We consider four representative scenarios and show the detailed problem setups and additional investigations as follows.

Manipulating mechanical stress distribution in varied geometric regions
The design domain is 180 × 90 mm 2 , and we use a structured background mesh consisting of 30 × 15 rectangular elements for Cases 1 and 2 and use an unstructured background mesh with a similar element size for Case 3. We set the maximum optimization iteration as 150 and p = 2 in the objective function.In the main manuscript, we modulate the hydrostatic stress to a target value of σ h = 0.28 MPa under an applied displacement of u = 1.5 mm.The selection of this setup is solely for demonstration purposes, and the proposed computational framework is not confined to a specific stress measure, target value, or applied displacement.To illustrate the versatility of the framework, we take Case 2 as an example and explore stress modulation effects under varying setups.Supplementary Figure 6 displays the initial, optimized, and target stresses, with the insets illustrating the design setup.Specifically, in Supplementary Figure 6a, we modulate the shear stress (instead of hydrostatic stress) to two target values, τ = 0.02 MPa and τ = 0.04 MPa, respectively.In Supplementary Figure 6b, we explore two different target values for hydrostatic stress, σ h = 0.18 MPa and σ h = 0.38 MPa (as opposed to σ h = 0.28 MPa).Finally, in Supplementary Figure 6c, we consider two different applied displacements, u = 0.20 mm and u = 0.25 mm (instead of u = 0.15 mm).Supplementary Table 1 summarizes the design setups and the corresponding relative stress modulation errors.According to Supplementary Figure 6 and Supplementary Table 1, the stress distribution aligns with the target after optimization across diverse setups, indicating the robust generality of the proposed framework for stress modulation.

Simultaneous mechanical stress modulation in multiple complex regions with distinct target values
The design domain is 350 × 100 mm 2 , and we use a structured background mesh consisting of 875 × 250 rectangular elements.We set the maximum optimization iteration as 200 and p = 2 in the objective function.

Bio-inspired virtual growth simulator
Based on the optimized physical density and frequency fields obtained from macroscopic topology optimization, we can utilize a bio-inspired virtual growth simulator -extended from the original virtual growth scheme [1] -to generate irregular architected materials with heterogeneous disordered microstructures.In this section, we introduce the basic procedures of this virtual growth simulator and the involved two-mesh-projection scheme.

Basic procedures of the virtual growth simulator
During the macroscopic topology optimization, we discretize a design domain with a (either structured or unstructured) background mesh and obtain the optimized density and frequency fields defined on this background mesh.Supplementary Figure 9a shows one representative background mesh consisting of four grids: one solid grid (ρ = 1), one void grid (ρ = 0), and two designable grids associated with an optimized frequency combination, {ξ 1 , ξ 2 , ..., ξ N b }.To utilize these optimized fields for guiding the material generation, we propose a virtual growth simulator consisting of the following basic procedures: defining a structured foreground mesh fully covering the background mesh (Supplementary Figure 9b), projecting the optimized density and frequency fields onto the defined foreground mesh (Supplementary Figure 9c; see details in Section 4.2), growing disordered microstructures based on the projected fields, and cutting the foreground mesh to fit the geometry of the background mesh (Supplementary Figure 9d).Finally, we obtain seamless integration of irregular architected materials with optimized heterogeneous disordered microstructures (Supplementary Figure 9d).

Two-mesh-projection scheme
As mentioned above, we need to project the optimized density and frequency fields onto a structured foreground mesh within the virtual growth simulator (Supplementary Figure 9b-c).This is because the original virtual growth scheme [1] is only applicable for a structured mesh with either square (in 2D) or cubic (in 3D) elements.To accomplish this goal, we propose a two-mesh-projection scheme to determine the corresponding source background grid for each foreground one.Specifically, we calculate the area (in 2D) or volume (in 3D) ratio between each pair of the intersected region (between background and foreground grids) and the foreground grid.The foreground grid then inherits the optimized density and frequency fields of the background grid associated with the highest area or volume ratio.In addition, a foreground grid is classified as solid if it is entirely within one or more solid background grids; a foreground grid is classified as void if it does not intersect with any designable background grids and is not classified as solid.These rules ensure that the obtained irregular architected materials are still fully connected after we cut them to match the geometries of the background mesh (Supplementary Figure 9d).
The pseudo-code of the proposed two-mesh-projection scheme is given in Algorithm 1.The proposed scheme involves computing the area ratios between polygons (in 2D) or volume ratios between polyhedra (in 3D) for each pair of the intersected region and the foreground grid.To accomplish this goal, we subdivide each background and foreground grid into many subgrids.We then calculate the central coordinates of each subgrid and use these coordinates to determine the nearest background subgrid for each foreground subgrid.By associating each foreground subgrid with the background grid containing the nearest background subgrid, we can approximate the area or volume ratios between the intersected region and the foreground grid.To effectively perform these complex geometric operations, we use an in-house Python code [5] with the PyVista [18], GeoPandas [19], and Shapely [20] packages.Notably, our implementation has no restrictions on the geometries of both the background and foreground meshes, enabling the projection of optimized fields from a structured/unstructured background mesh onto another structured/unstructured foreground mesh.
Algorithm 1: Two-mesh-projection scheme for projecting the optimized fields 1 Inputs: the number of foreground grids, N fore ; the number of background subgrids in one background grid, N back sub ; the number of foreground subgrids in one foreground grid, N fore sub ; 2 Find all solid and void background grids, S and V, respectively; 3 Subdivide each background grid into N back sub subgrids; 4 Subdivide each foreground grid into N fore sub subgrids; 5 Initialize the foreground grid count, n fore = 1; 6 while n fore ≤ N fore do 7 Initialize a relation vector, R = 0 ∈ R N fore sub .Here R i = j represents that subgrid i of foreground grid n fore is within background grid j; subgrid i does not belong to any background grid if R i = 0;

Average stress measures
To conduct a stress convergence study for irregular architected materials, we need to determine an appropriate measure for evaluating the average stress.We start from the stress tensor, σ(X), averaged over background grid e, which can be expressed as [21]

(X)dX
where σ e is the average stress tensor, and • is the operator for computing the volume average.The variable Ω e is the domain occupied by background grid e.We now consider an arbitrary stress measure, σ s = f (σ), as a function of the stress tensor, σ.The average stress measure can be computed as and it holds σ s e = f (σ) e = σ s e if the stress measure, σ s , is a linear combination of the components (σ ij for i, j = 1, 2, 3) of the stress tensor, σ.
In this study, we mainly focus on the hydrostatic stress and the shear stress, and both of them are linear functions of the stress components.Therefore, the macroscopic average stress measure (Supplementary Figure 9a) is given as σ s,macro e = σ s,macro e = 1 |Ω e | Ωe σ s,macro (X)dX, (5) and the microscopic average stress measure (Supplementary Figure 9d) is given as where σ s,macro and σ s,micro are macroscopic and microscopic stress measures of interest, respectively.

Supplementary Figure 2 :
Material property distribution of the created specimens.The x-axis is the material property (D11 component of the D matrix here), the left y-axis is the statistical frequency, and right y-axis is the probability density function.The inset color wheels represent the frequency combinations of building blocks.

Supplementary Figure 4 : 5 :
Architecture of a fully connected neural network for predicting homogenized elastic modulus.The architecture contains one input layer with N b nodes corresponding to the N b entries of the input frequency combination.It also contains two hidden layers with m and p nodes, respectively, and one output layer with 6 nodes in 2D and 21 nodes in 3D corresponding to the maximum number of independent entries of homogenized elastic modulus.Determining the optimal architecture parameters and formally training the neural network.a, Using a stochastic search to determine the optimal architecture parameters of the neural network.This figure shows the histories of the loss functions of 100 sets of randomly generated architecture parameters, including the learning rate of the optimizer and the numbers of nodes of the two hidden layers.b, The complete history of the loss function during the formal training of weighting matrices and bias vectors.

Supplementary Figure 6 :
Investigation of stress modulation effects under varying setups.a, The stress measure is changed from hydrostatic to shear stress.b, The target value is changed from 0.28 MPa to 0.18 and 0.38 MPa.c, The applied displacement is changed from 0.15 mm to 0.20 and 0.25 mm.The insets illustrate the updated design setups.

8 :
Finite element discretization and the optimized displacement field in z-direction of orthopedic femur restoration.The model consists of the femur embedded with a pre-fracture and artificial support.Both components are discretized with tetrahedral elements due to complex geometries.In the visualization, color signifies the optimized displacement field in the z-direction, with the inset offering a zoom-in view of the displacement field around the fracture.

Supplementary Figure 9 :
The proposed virtual growth simulator.a, One representative unstructured background mesh consisting of four elements.b, A structured foreground mesh defined to cover the unstructured background mesh.c, The structured foreground mesh with projected density and frequency fields.d, Optimized irregular architected materials cut to match the geometries of the background mesh.

21 else 22
Remove all entries in R if they are also in S ∪ V ∪ {0}; 23 Foreground grid n fore inherits the optimized density and frequency fields of background grid n = Mode(R);

5. 2 .Supplementary Figure 10 :
Convergence of the average stress b a Stress convergence studies for mechanical stress modulation within distinct control regions.a, Case 1; b, Case 2. The plots illustrate the relationship between the hydrostatic stress measure, denoted by σ h (in MPa), and the number of basic building blocks in one direction within one microstructure.Each error bar represents one numerical sample.The circular dot of the error bar represents the mean hydrostatic stress in the control region, and the half length of the error bar represents the standard deviation.
h 1,m ] , and h 2 = [h 2,1 , h 2,2 , ..., h 2,p ] , are intermediate vectors used in the forward pass of the neural network.Weighting matrices, W 1 , W 2 , and W 3 , and bias vectors, b 1 , b 2 , and b 3 , are parameters to be optimized such that the prediction error of microstructural elastic modulus, D vec (the vector form of D), is minimized.Superscripts indicate the matrix and vector dimensions for the reader's convenience.
Supplementary Table 1: Summary of design setups and relative stress errors in Supplementary Figure6

8
Initialize the foreground subgrid count, n fore sub = 1;Find the nearest background subgrid and the associated background grid, k; 15 end 16 Let U denote unique entries in R; 17 if U ⊆ S then 18 Foreground grid n fore is solid; 19 else if U ⊆ S ∪ V ∪ {0} then 20 Foreground grid n fore is void;